#########################################################
#########################################################
####     R Code for German Colonialism Project      #####
####       Goal: Create Major Figures               #####
####                Jan Pierskalla                  #####
####                  09/06/2016                    ##### 
#########################################################
#########################################################

###########################################
#### 1. Load Packages and Functions
###########################################

rm(list=ls(all=TRUE))

library(raster)
library(maptools)
library(rgdal)
library(raster)
library(rgeos)
library(lattice)
library(fields)
library(gplots)
library(spdep)
library(cshapes)
library(countrycode)
library(sp)
library(rgdal)
library(foreign)
library(ggplot2)
library(RColorBrewer)
library(animation)

plot.heat <- function(counties.map,z,title=NULL,breaks=NULL,reverse=FALSE,cex.legend=1,bw=.2,col.vec=NULL,plot.legend=TRUE) {
  ##Break down the value variable
  if (is.null(breaks)) {
    breaks=
      seq(
        floor(min(counties.map@data[,z],na.rm=TRUE)*10)/10
        ,
        ceiling(max(counties.map@data[,z],na.rm=TRUE)*10)/10
        ,.1)
  }
  counties.map@data$zCat <- cut(counties.map@data[,z],breaks,include.lowest=TRUE)
  cutpoints <- levels(counties.map@data$zCat)
  if (is.null(col.vec)) col.vec <- gray.colors(length(levels(counties.map@data$zCat)))
  if (reverse) {
    cutpointsColors <- rev(col.vec)
  } else {
    cutpointsColors <- col.vec
  }
  levels(counties.map@data$zCat) <- cutpointsColors
  plot(counties.map,border="transparent",lwd=0.00001,axes = FALSE, las = 1,col=as.character(counties.map@data$zCat))
  #if (!is.null(state.map)) {
  # plot(state.map,add=TRUE,lwd=1)
  #}
  ##with(counties.map.c,text(x,y,name,cex=0.75))
  if (plot.legend) legend("bottomleft", cutpoints, fill = cutpointsColors,bty="n",title=title,cex=cex.legend)
  ##title("Cartogram")
}


###########################################
#### 2. Set Working Directory
###########################################


setwd("~/Dropbox/German_Colonialism/Statebuilding_Paper/Replication_Archive/")



###########################################
#### 3. Make Maps
###########################################


# Plot of stations, by year
load("~/Dropbox/German_Colonialism/Statebuilding_Paper/Replication_Archive/1890_grid_stations.Rdata")
grid_stations@data$id = rownames(grid_stations@data)
grid_stations_df = fortify(grid_stations,region="id")
grid_stations_df = join(grid_stations_df, grid_stations@data, by="id")
grid_stations_df$Stations <- ifelse(grid_stations_df$station_count_nomiss==1,"Station","No Station")
gs.pal <- colorRampPalette(c("#f0f0f0","#636363"),bias=.1,space="rgb")
pdf(file="~/Dropbox/German_Colonialism/Statebuilding_Paper/Figures/Stations_1890.pdf",width=10)
map <- qplot(long, lat, data=grid_stations_df, group=group , fill=Stations, geom="polygon")+
  scale_fill_manual(values=gs.pal(2))+ coord_equal()+geom_path(size=0.1,color="black") +theme_bw()+ ggtitle("Colonial Stations, 1890")
map
dev.off()

load("~/Dropbox/German_Colonialism/Statebuilding_Paper/Replication_Archive/1909_grid_stations.Rdata")
grid_stations@data$id = rownames(grid_stations@data)
grid_stations_df = fortify(grid_stations,region="id")
grid_stations_df = join(grid_stations_df, grid_stations@data, by="id")
grid_stations_df$Stations <- ifelse(grid_stations_df$station_count_nomiss==1,"Station","No Station")
gs.pal <- colorRampPalette(c("#f0f0f0","#636363"),bias=.1,space="rgb")
pdf(file="~/Dropbox/German_Colonialism/Statebuilding_Paper/Figures/Stations_1909.pdf",width=10)
map <- qplot(long, lat, data=grid_stations_df, group=group , fill=Stations, geom="polygon")+
  scale_fill_manual(values=gs.pal(2))+ coord_equal()+geom_path(size=0.1,color="black") +theme_bw()+ ggtitle("Colonial Stations, 1909")
map
dev.off()



##############################################
#### Simulation of Substantive Effect
##############################################
library(ggplot2)
library(gpclib)
library(maptools)
library(gdata)
library(Zelig)
library(scales)
library(foreign)
library(plyr)
library(WDI)
library(cshapes)
library(sp)
library(classInt)
library(reshape2)
library(MASS)
library(lmtest)

#Battle_Index Effect
data <- read.dta("~/Dropbox/German_Colonialism/Statebuilding_Paper/Replication_Archive/grid_panel_replication.dta")
effect <- read.dta("~/Dropbox/German_Colonialism/Statebuilding_Paper/Replication_Archive/effect_battle.dta")

datavec <- seq(-2.3,2.2,0.1)
plots <- data.frame(datavec,effect$estimate)
plotdata <- data.frame(data$battle_index_l)
colnames(plotdata)[1] <- "battle_index_l"

pdf("~/Dropbox/German_Colonialism/Statebuilding_Paper/Replication_Archive/Figures/effect_battle.pdf",height=7, width=7*1.62)
ggplot(plots,aes(datavec,effect$estimate))+geom_line()+
  theme_bw()+scale_y_continuous(limits = c(0,0.25))+
  theme(axis.text=element_text(size=17),axis.title=element_text(size=17,face="bold"))+
  xlab("Battle Index")+ylab("P(Station)")+geom_ribbon(aes(ymin=effect$min95, ymax=effect$max95),alpha=0.2)+
  geom_point(data=plotdata, aes(x=plotdata$battle_index_l, y=0), alpha=0.2)
dev.off()  



# stratval_n_correct_l

data <- read.dta("~/Dropbox/German_Colonialism/Statebuilding_Paper/Replication_Archive/grid_panel_replication.dta")
effect <- read.dta("~/Dropbox/German_Colonialism/Statebuilding_Paper/Replication_Archive/effect_strat.dta")

datavec <- seq(-0.655,5.3419,0.1)
plots <- data.frame(datavec,effect$estimate)
plotdata <- data.frame(data$stratval_n_correct_l)
colnames(plotdata)[1] <- "stratval_n_correct_l"

pdf("~/Dropbox/German_Colonialism/Statebuilding_Paper/Replication_Archive/Figures/effect_strat.pdf",height=7, width=7*1.62)
ggplot(plots,aes(datavec,effect$estimate))+geom_line()+
  theme_bw()+scale_y_continuous(limits = c(0,1))+
  theme(axis.text=element_text(size=17),axis.title=element_text(size=17,face="bold"))+
  xlab("Territorial Control Value")+ylab("P(Station)")+geom_ribbon(aes(ymin=effect$min95, ymax=effect$max95),alpha=0.2)+
  geom_point(data=plotdata, aes(x=plotdata$stratval_n_correct_l, y=0), alpha=0.2)
dev.off()  




########################################################
########################################################
# Create Figure for Territorial Value Variable
########################################################
########################################################

load("~/Dropbox/German_Colonialism/Statebuilding_Paper/Replication_Archive/1890_grid_stations.Rdata")

plot(grid_stations)
grid_stations@data$station_count[is.na(grid_stations@data$station_count)] <- 0

stations <- coordinates(grid_stations[grid_stations@data$station_count==1,])
grid_centroids <- coordinates(grid_stations)
distances_station <- spDists(grid_centroids,stations,longlat=TRUE)
distances_station <- data.frame(distances_station)
rownames(distances_station) <- grid_stations@data$pu_id
colnames(distances_station) <- grid_stations@data$pu_id[grid_stations@data$station_count==1]

min <- NA
mindist_id <- NA

for (i in 1:681){
  distances <- c(t(distances_station[i,]))
  ids <- c(as.numeric(colnames(distances_station[i,])))
  mindist <- data.frame(distances, ids)
  min[i] <- min(mindist$distances,na.rm=T)
  mindist_id[i] <- mindist$ids[mindist==min[i]]
}

dfPoints <- matrix(NA,nrow=681,ncol=4)
dfPoints <- data.frame(dfPoints)
colnames(dfPoints) <- c("xlat","xlon","ylat","ylon")

stations <- coordinates(grid_stations[grid_stations@data$station_count==1,])
colnames(stations) <- c("lon","lat")
stations <- data.frame(stations)
stations$id <- grid_stations@data$pu_id[grid_stations@data$station_count==1]

for (i in 1:681){
  
  dfPoints$xlon[i] <- grid_centroids[i,1]
  dfPoints$xlat[i] <- grid_centroids[i,2]
  
  dfPoints$ylon[i] <- stations$lon[stations$id==mindist_id[i]]
  dfPoints$ylat[i] <- stations$lat[stations$id==mindist_id[i]]
}

pdf(file="~/Dropbox/German_Colonialism/Statebuilding_Paper/Replication_Archive/Figures/measure1b.pdf",width=10)
dist <- c(-1,0,1)
plot.heat(grid_stations,z="station_count",breaks=dist,reverse=TRUE,plot.legend=F)
segments(x0=dfPoints$xlon,y0=dfPoints$xlat,x1=dfPoints$ylon,y1=dfPoints$ylat,col="darkblue")
dev.off()

pdf(file="~/Dropbox/German_Colonialism/Statebuilding_Paper/Replication_Archive/Figures/measure1a.pdf",width=10)
dist <- c(-1,0,1)
plot.heat(grid_stations,z="station_count",breaks=dist,reverse=TRUE,plot.legend=F)
#segments(x0=dfPoints$xlon,y0=dfPoints$xlat,x1=dfPoints$ylon,y1=dfPoints$ylat,col="darkblue")
dev.off()




# now with a hypothetical new station
load("~/Dropbox/German_Colonialism/Statebuilding_Paper/Replication_Archive/1890_grid_stations.Rdata")

plot(grid_stations)
grid_stations@data$station_count[is.na(grid_stations@data$station_count)] <- 0

grid_stations@data$station_count[530] <- 1

dist <- c(-1,0,1)
plot.heat(grid_stations,z="station_count",breaks=dist,reverse=TRUE,plot.legend=F)

stations <- coordinates(grid_stations[grid_stations@data$station_count==1,])
grid_centroids <- coordinates(grid_stations)
distances_station <- spDists(grid_centroids,stations,longlat=TRUE)
distances_station <- data.frame(distances_station)
rownames(distances_station) <- grid_stations@data$pu_id
colnames(distances_station) <- grid_stations@data$pu_id[grid_stations@data$station_count==1]

min <- NA
mindist_id <- NA

for (i in 1:681){
  distances <- c(t(distances_station[i,]))
  ids <- c(as.numeric(colnames(distances_station[i,])))
  mindist <- data.frame(distances, ids)
  min[i] <- min(mindist$distances,na.rm=T)
  mindist_id[i] <- mindist$ids[mindist==min[i]]
}

dfPoints <- matrix(NA,nrow=681,ncol=4)
dfPoints <- data.frame(dfPoints)
colnames(dfPoints) <- c("xlat","xlon","ylat","ylon")

stations <- coordinates(grid_stations[grid_stations@data$station_count==1,])
colnames(stations) <- c("lon","lat")
stations <- data.frame(stations)
stations$id <- grid_stations@data$pu_id[grid_stations@data$station_count==1]

for (i in 1:681){
  
  dfPoints$xlon[i] <- grid_centroids[i,1]
  dfPoints$xlat[i] <- grid_centroids[i,2]
  
  dfPoints$ylon[i] <- stations$lon[stations$id==mindist_id[i]]
  dfPoints$ylat[i] <- stations$lat[stations$id==mindist_id[i]]
}

pdf(file="~/Dropbox/German_Colonialism/Statebuilding_Paper/Replication_Archive/Figures/measure2b.pdf",width=10)
dist <- c(-1,0,1)
plot.heat(grid_stations,z="station_count",breaks=dist,reverse=TRUE,plot.legend=F)
segments(x0=dfPoints$xlon,y0=dfPoints$xlat,x1=dfPoints$ylon,y1=dfPoints$ylat,col="darkblue")
dev.off()

myColours <- rep(NA, 681)
myColours[530] <- "red"

pdf(file="~/Dropbox/German_Colonialism/Statebuilding_Paper/Replication_Archive/Figures/measure2a.pdf",width=10)
dist <- c(-1,0,1)
plot.heat(grid_stations,z="station_count",breaks=dist,reverse=TRUE,plot.legend=F)
plot(grid_stations,col=myColours,add=T,border="transparent",lwd=0.00001)
#segments(x0=dfPoints$xlon,y0=dfPoints$xlat,x1=dfPoints$ylon,y1=dfPoints$ylat,col="darkblue")
dev.off()

